So, to start out I need to get some data. I’m thinking of taking; - California Counties data (X) * US boundaries data package has it. - average rainfall (X) * got it from national weather sevice. - county population(X) * got it from census - fire polygons or rasters (X) * I found exactly what I wanted. - city names (X) * I ended up getting this from another source. - city polygons (X) * Got it from the city names source.

and basically combining all of these into an interractive map that the user can zoom in using the leaflet package, flip through the years and also see a graph showing the estimated number of people whose lives were affected by the fires. I could also add other fitting data if I find good sources, but realistically, I might struggle to even find the data that I have listed here.

Tomorrow I’ll get started tracking down the data.

library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ ggplot2 3.3.4     ✔ purrr   0.3.4
## ✔ tibble  3.1.2     ✔ dplyr   1.0.7
## ✔ tidyr   1.1.3     ✔ stringr 1.5.0
## ✔ readr   1.4.0     ✔ forcats 0.5.1
## Warning: package 'stringr' was built under R version 4.1.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
library(sf)
## Warning: package 'sf' was built under R version 4.1.2
## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(rgdal)
## Warning: package 'rgdal' was built under R version 4.1.2
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-3, (SVN revision 1187)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.4.2, released 2022/03/08
## Path to GDAL shared files: /Users/Jason/Library/R/x86_64/4.1/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: /Users/Jason/Library/R/x86_64/4.1/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
GDALinfo("../data/nws_precip_ytd_20220101_conus.tif")
## Warning: GDAL support is provided by the sf and terra packages among others
## Warning: GDAL support is provided by the sf and terra packages among others
## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver

## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver

## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver

## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver
## Warning: GDAL support is provided by the sf and terra packages among others
## rows        881 
## columns     1121 
## bands       4 
## lower left origin.x        -1904912 
## lower left origin.y        -7619987 
## res.x       4763 
## res.y       4763 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-105 +x_0=0 +y_0=0 +R=6371200
## +units=m +no_defs +type=crs 
## file        ../data/nws_precip_ytd_20220101_conus.tif 
## apparent band summary:
##    GDType hasNoDataValue   NoDataValue blockSize1 blockSize2
## 1 Float32           TRUE -3.402823e+38        256        256
## 2 Float32           TRUE -3.402823e+38        256        256
## 3 Float32           TRUE -3.402823e+38        256        256
## 4 Float32           TRUE -3.402823e+38        256        256
## apparent band statistics:
##          Bmin       Bmax Bmean Bsd
## 1 -4294967295 4294967295    NA  NA
## 2 -4294967295 4294967295    NA  NA
## 3 -4294967295 4294967295    NA  NA
## 4 -4294967295 4294967295    NA  NA
## Metadata:
## AREA_OR_POINT=Area 
## creation_time=20220112T010818 
## data_time=20220101
precipitation_conus_path <- "../data/nws_precip_ytd_20220101_conus.tif"

precip_conus_test <- raster(precipitation_conus_path, band = 2)

spplot(precip_conus_test)

Okay, so I’ve tentatively got my precipitation data, wait, let’s check that.

Here’s the precipitation data source

https://water.weather.gov/precip/download.php

Source: https://www.noaa.gov/

GeoTIFF

The new QPE GeoTIFFs generated from the NCEP Stage IV data are multi-band GeoTIFF. The bands they contain are:

Band 1 - Observation - Last 24 hours of QPE spanning 12Z to 12Z in inches
Band 2 - PRISM normals - PRISM normals in inches (see "Normal Precipitation" section on the About page)
Band 3 - Departure from normal - The departure from normal in inches
Band 4 - Percent of normal - The percent of normal

Now, this raster shows precip in inches of rainfall. I’m not going to be using this for a quantitative calculation, so simply having this much is enough for me.

The next step will be downloading all of the yearly data rasters from 2016-2022, then, once I actually have all of them, I have to figure out how to isolate the California part of the data. Most likely, I will be projecting the state boundary onto the raster and then selecting the raster points within the polygon. Should be simple enough once I actually get around to it. I need to also find the rest of my data.

The data protocol changed in 2017 for the precip data, so I might need to do something different for year 2016 precip data.

Now to check on the r package data sources.

library(USAboundaries)

california_state_select <- us_boundaries(states = 'california')
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
cali_counties_data <- us_boundaries(states = 'california', type = "county")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
plot(california_state_select)
## Warning: plotting the first 9 out of 12 attributes; use max.plot = 12 to plot
## all

plot(cali_counties_data)
## Warning: plotting the first 9 out of 12 attributes; use max.plot = 12 to plot
## all

Okay, so I’ve got the cali state and county polys in this package. Next is to check on the uscities.csv, and find out where its data comes from as well as the us boundary data comes from.

Actually, there is no data source imbeded in the csv that I can find so I’m gonna try some open source data.

https://data.ca.gov/dataset/ca-geographic-boundaries

Source: https://data.ca.gov/

this is the one I’ll try

library(sf)

ca_places_path <- "../data/ca-places-boundaries/CA_Places_TIGER2016.shp"

ca_places_bound_test <- st_read(ca_places_path) %>% 
  st_as_sf() %>% 
  filter(NAME == "Santa Rosa")
## Reading layer `CA_Places_TIGER2016' from data source 
##   `/Users/Jason/github/fire_map_proj/data/ca-places-boundaries/CA_Places_TIGER2016.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1522 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13833620 ymin: 3833636 xmax: -12715920 ymax: 5159958
## Projected CRS: WGS 84 / Pseudo-Mercator
plot(ca_places_bound_test)
## Warning: plotting the first 9 out of 16 attributes; use max.plot = 16 to plot
## all

Okay, this looks promising. I’m just not 100% sure if it will give me city polygons or what it even is really. It actually looks like exactly what I wanted.

Now onto fire data for california. I’ll try this dataset first.

https://hub-calfire-forestry.hub.arcgis.com/datasets/CALFIRE-Forestry::fire-perimeters/explore?location=37.417017%2C-122.003310%2C10.00

Source: https://www.fire.ca.gov/

Seems promising. Time to start unpacking the data.

fire_data_path <- "../data/California_Wildland_Fire_Perimeters_-All/California_Wildland_Fire_Perimeters_All.shp"

fire_data_perim_test <- st_read(fire_data_path) %>% 
  st_as_sf() %>% 
  filter(FIRE_NAME == "TUBBS")
## Reading layer `California_Wildland_Fire_Perimeters_All' from data source 
##   `/Users/Jason/github/fire_map_proj/data/California_Wildland_Fire_Perimeters_-All/California_Wildland_Fire_Perimeters_All.shp' 
##   using driver `ESRI Shapefile'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received an unexpected geometry. Either a polygon
## with interior rings, or a polygon with less than 4 points, or a non-Polygon
## geometry. Return arguments as a collection.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Geometry of polygon of fid 20145 cannot be translated to Simple
## Geometry. All polygons will be contained in a multipolygon.
## Simple feature collection with 21688 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13848330 ymin: 3833204 xmax: -12705610 ymax: 5255380
## Projected CRS: WGS 84 / Pseudo-Mercator
print(fire_data_perim_test)
## Simple feature collection with 1 feature and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13666000 ymin: 4645848 xmax: -13646670 ymax: 4674086
## Projected CRS: WGS 84 / Pseudo-Mercator
##   OBJECTID YEAR_ STATE AGENCY UNIT_ID FIRE_NAME  INC_NUM ALARM_DATE  CONT_DATE
## 1    41660  2017    CA    CDF     LNU     TUBBS 00010045 2017-10-08 2017-11-01
##   CAUSE COMMENTS REPORT_AC GIS_ACRES C_METHOD OBJECTIVE FIRE_NUM Shape__Are
## 1    14     <NA>     36701  36701.98        3         1     <NA>  243263474
##   Shape__Len          COMPLEX_NA COMPLEX_IN                       geometry
## 1   167697.4 CENTRAL LNU COMPLEX   00010104 MULTIPOLYGON (((-13649169 4...
plot(fire_data_perim_test)
## Warning: plotting the first 10 out of 20 attributes; use max.plot = 20 to plot
## all

Okay, that was great. Now I need to get the populations of each county. I found a data source here. It’s xslx format, but that should be fine.

https://www.census.gov/data/datasets/time-series/demo/popest/2010s-counties-total.html

Source: https://www.census.gov/en.html

library(readxl)

county_pop_path <- "../data/co-est2019-annres-06.xlsx"

county_pop <- read_excel(county_pop_path)
## New names:
## • `` -> `...2`
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`

Okay, looks good mostly. The csv doesn’t have good column names, but that’s fine. They’re still labeled, just not conveniently so I can work with it. So, these numbers are the estimates from between 2010 and 2019 which is basically just a rough estimation which is close enough for me. These estimates are more than good enough for what I intend to use them for.

OKAY. So I now have all my data. What is the next step? I should get started with making one of these data sets my ‘base’ and adding the data layers on from there. What is my core data layer? …

It’s probably best to start with the US boundaries data set. It will give me the state and county polygons. I also need to figure out what data is getting joined to what. So, there are going to be a few separate layers. The base layer will be the leaflet map that I’m projecting the data onto. Just above the base layer should be the precipitation layer so as to keep the polygon layer boundaries more clearly defined. The layer above that will be the county layer which will have the county population data joined to it. Then, the next layer up will be the city polygon layer and the fire boundary layer on top. These two layers will help illustrate the proximity of the fires to major city centers.

With the map breakdown like this, I can create the graph of population ‘affected’ by fires relatively easily.

So in order to get everything ready for the map, I need to start cleaning up the data and joining what needs to be joined. To do this, I should make a ggplot that allows me to visualize the data layers from bottom to top and make sure I’m cutting out any data entries or columns that I don’t need.

Let’s get started with the us boundaries ggplot

#libraries are already called, but I might need to reimport the bounradies

#get the california state boundary
california_state <- us_boundaries(states = 'california')
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
#get the california counties as an sfc
cali_counties <- us_boundaries(states = 'california', type = "county")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
# I didn't need to recall the boundaries, but I'll just leave them here for now.

#okay, now I want to join the population data to cali_counties.
#we want a left join where pop gets added onto the states based on name.
print(county_pop)
## # A tibble: 67 x 13
##    `table with row h… ...2   ...3  ...4     ...5    ...6    ...7    ...8    ...9
##    <chr>              <chr>  <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 Annual Estimates … <NA>   <NA>  <NA>  NA      NA      NA      NA      NA     
##  2 Geographic Area    40269  <NA>  Popu… NA      NA      NA      NA      NA     
##  3 <NA>               Census Esti… 2010   2.01e3  2.01e3  2.01e3  2.01e3  2.02e3
##  4 California         37253… 3725… 3731…  3.76e7  3.79e7  3.83e7  3.86e7  3.89e7
##  5 .Alameda County, … 15102… 1510… 1512…  1.53e6  1.55e6  1.58e6  1.61e6  1.63e6
##  6 .Alpine County, C… 1175   1175  1161   1.09e3  1.11e3  1.13e3  1.08e3  1.08e3
##  7 .Amador County, C… 38091  38091 37886  3.75e4  3.71e4  3.66e4  3.67e4  3.70e4
##  8 .Butte County, Ca… 220000 2200… 2199…  2.20e5  2.21e5  2.22e5  2.24e5  2.25e5
##  9 .Calaveras County… 45578  45578 45468  4.52e4  4.48e4  4.47e4  4.47e4  4.50e4
## 10 .Colusa County, C… 21419  21407 21437  2.13e4  2.13e4  2.12e4  2.12e4  2.12e4
## # … with 57 more rows, and 4 more variables: ...10 <dbl>, ...11 <dbl>,
## #   ...12 <dbl>, ...13 <dbl>
print(cali_counties$name)
##  [1] "San Francisco"   "Nevada"          "Merced"          "Alpine"         
##  [5] "Contra Costa"    "Tulare"          "Lake"            "Mono"           
##  [9] "Imperial"        "Alameda"         "Santa Clara"     "Stanislaus"     
## [13] "Lassen"          "Santa Cruz"      "Sierra"          "Napa"           
## [17] "Colusa"          "Sacramento"      "Sutter"          "Mariposa"       
## [21] "Humboldt"        "Glenn"           "Amador"          "Madera"         
## [25] "Tuolumne"        "Shasta"          "Plumas"          "Yolo"           
## [29] "Kern"            "Fresno"          "Kings"           "San Joaquin"    
## [33] "Los Angeles"     "Ventura"         "Santa Barbara"   "Marin"          
## [37] "San Bernardino"  "Inyo"            "Modoc"           "Yuba"           
## [41] "San Benito"      "Tehama"          "Siskiyou"        "Del Norte"      
## [45] "Orange"          "Placer"          "El Dorado"       "San Diego"      
## [49] "Calaveras"       "Monterey"        "Riverside"       "Butte"          
## [53] "San Mateo"       "Sonoma"          "Trinity"         "Mendocino"      
## [57] "Solano"          "San Luis Obispo"
#so, it appears that the names don't match. Trying to remove the .~~~ County, California
#could be annoying, so I might just add that onto the cali_counties name and join by that
#I can try mutating on a new column with altered nameing.
cali_counties_join <- cali_counties %>% 
  mutate(name_join = paste( ".", cali_counties$name, " County, California", sep = ""))

print(cali_counties_join)
## Simple feature collection with 58 features and 13 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -124.4096 ymin: 32.53416 xmax: -114.1391 ymax: 42.00925
## Geodetic CRS:  WGS 84
## First 10 features:
##     statefp countyfp countyns       affgeoid geoid          name lsad
## 76       06      075 00277302 0500000US06075 06075 San Francisco   06
## 371      06      057 01682927 0500000US06057 06057        Nevada   06
## 636      06      047 00277288 0500000US06047 06047        Merced   06
## 714      06      003 01675840 0500000US06003 06003        Alpine   06
## 717      06      013 01675903 0500000US06013 06013  Contra Costa   06
## 774      06      107 00277318 0500000US06107 06107        Tulare   06
## 803      06      033 00277281 0500000US06033 06033          Lake   06
## 898      06      051 00277290 0500000US06051 06051          Mono   06
## 953      06      025 00277277 0500000US06025 06025      Imperial   06
## 982      06      001 01675839 0500000US06001 06001       Alameda   06
##           aland    awater state_name state_abbr jurisdiction_type
## 76    121455687 479136515 California         CA             state
## 371  2480606160  41513172 California         CA             state
## 636  5012175306 112509475 California         CA             state
## 714  1912292633  12557304 California         CA             state
## 717  1871930816 209819213 California         CA             state
## 774 12495016292  37081410 California         CA             state
## 803  3254288286 188912510 California         CA             state
## 898  7896839510 214694734 California         CA             state
## 953 10817291640 790230304 California         CA             state
## 982  1914242789 212979931 California         CA             state
##                           geometry                         name_join
## 76  MULTIPOLYGON (((-122.512 37... .San Francisco County, California
## 371 MULTIPOLYGON (((-121.2661 3...        .Nevada County, California
## 636 MULTIPOLYGON (((-121.2268 3...        .Merced County, California
## 714 MULTIPOLYGON (((-120.0725 3...        .Alpine County, California
## 717 MULTIPOLYGON (((-122.4253 3...  .Contra Costa County, California
## 774 MULTIPOLYGON (((-119.5461 3...        .Tulare County, California
## 803 MULTIPOLYGON (((-123.0778 3...          .Lake County, California
## 898 MULTIPOLYGON (((-119.6392 3...          .Mono County, California
## 953 MULTIPOLYGON (((-116.1056 3...      .Imperial County, California
## 982 MULTIPOLYGON (((-122.3337 3...       .Alameda County, California
#now try a left join, horrific column name
cali_counties_join <- cali_counties_join %>% 
  left_join(county_pop, by = c("name_join" = "table with row headers in column A and column headers in rows 3 through 4 (leading dots indicate sub-parts)"))

#great, that actually worked, now we can cut cali_counties down to size with select and filter.
#XXXXXXXXXX Continue from here XXXXXXXXXX
cali_counties_final <- cali_counties_join %>% 
  st_as_sf(crs = 4326) %>% 
  dplyr::select(name, '...2', geometry) %>%  #I seriously have to specify this? Completely ridiculous
  rename("county_population" = "...2") %>% 
  rename("county_name" = "name")

#create ggplot
ggplot()+
  #color the cali counties and their borders
  geom_sf(data = cali_counties_final, fill = "#454742", col = "#8d8c8b", size = .5)+
  #create the theme for the background. I'm reusing an old theme
  theme(panel.background = element_rect(fill = "#272727",
                                        colour = "#272727",
                                        size = .5,
                                        linetype = "dashed"),
        panel.grid.major = element_line(size = .5,
                                        color = "#3d3d3d",
                                        linetype = "solid"))+
  #title and subtitle. I might need another legend position for credits.
  theme(legend.position = "hide")+
    labs(title = "TITLE",
         subtitle = "SUBTITLE",
         x = "",
         y ="")

Okay, I managed to solve the github push issue by creating a new repo and making sure to .gitignore the data folder. Now I have a tidy looking county pop poly var.

Now I need to decide what to do next. I feel like I want to set up the leaflet now since the rest of the data just needs to be wrangled. Then again, getting the data tidy first would probably be for the best. Okay, I’ll focus on keeping the data tidy first and not get distracted.

#time to tidy the city polys
ca_places_path <- "../data/ca-places-boundaries/CA_Places_TIGER2016.shp"

ca_places_bound <- st_read(ca_places_path) %>% 
  st_as_sf(crs = 4326) %>% 
  dplyr::select(NAMELSAD, geometry) %>% 
  st_simplify(dTolerance = 200) %>%  #freeport ends up with 'polygon empty' keep note in case it becomes a problem later
  st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet
## Reading layer `CA_Places_TIGER2016' from data source 
##   `/Users/Jason/github/fire_map_proj/data/ca-places-boundaries/CA_Places_TIGER2016.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1522 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13833620 ymin: 3833636 xmax: -12715920 ymax: 5159958
## Projected CRS: WGS 84 / Pseudo-Mercator
#create ggplot
ggplot()+
  #color the objects
  geom_sf(data = ca_places_bound, fill = "#454742", col = "#8d8c8b", size = .5)+
  #create the theme for the background. I'm reusing an old theme
  theme(panel.background = element_rect(fill = "#272727",
                                        colour = "#272727",
                                        size = .5,
                                        linetype = "dashed"),
        panel.grid.major = element_line(size = .5,
                                        color = "#3d3d3d",
                                        linetype = "solid"))+
  #title and subtitle. I might need another legend position for credits.
  theme(legend.position = "hide")+
    labs(title = "TITLE",
         subtitle = "SUBTITLE",
         x = "",
         y ="")

Good, the next step should be to tidy up the fire map data and split it into different years. There’s probably a built in way of dealing with that, perhaps with groups or something, but I’m going to split the data for the readability of my code for now. I can pack things up nicer on the final map if I want to.

Okay, I seriously need to diagnose this terrible error. It’s so bad, that even the built in st_is_valid() function is useless and breaks upon this mighty error.

#fire map error diagnosis
fire_data_path <- "../data/California_Wildland_Fire_Perimeters_-All/California_Wildland_Fire_Perimeters_All.shp"

#okay, so the error I need to diagnose is that some of the fire shapefile geometries
#are broken and also completely break the sf::st_is_valid() diagnosis tool. I hate
#this and it's driving me insane, so I need to figure out how to isolate the issue
#without any actual tools. Let's start by recreating the error

#corrupt_fire_data_1 <- st_read(fire_data_path) %>%  #okay, st_read seems to work just fine
  #dplyr::select(YEAR_, FIRE_NAME, geometry) %>% #select also seems to work and shouldn't mess up anything
  #st_is_valid() %>% 
  #as.data.frame()
  #so, if I do the valid check now, the problem is still the same, all NA
  #what could I try next that might help isolate the issue? I can repeat the valid
  #check for all years I want to see if it's just 2017
  #filter(YEAR_ == 2016) %>% 
  #st_is_valid() %>% 
  #as.data.frame()

  #okay, so 2016 is still fine and works, let's try 2017
  #filter(YEAR_ == 2017) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #Once again, all NA

  #now try 2018
  #filter(YEAR_ == 2018) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #2018 is fine

  #filter(YEAR_ == 2019) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #2019 works too

  #filter(YEAR_ == 2020) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #2020 is fine

  #filter(YEAR_ == 2021) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #2021 isn't broken either

  #okay,so I have now established that the error is in 2017, now I'll create a new
  #chunk to keep my notes intact.

So, new chunk is going to be for testing the different ways I can break down the data from 2017 to locate the erroneous polygon that is breaking st_is_valid()

fire_data_path <- "../data/California_Wildland_Fire_Perimeters_-All/California_Wildland_Fire_Perimeters_All.shp"


corrupt_fire_data_1 <- st_read(fire_data_path) %>%  #okay, st_read seems to work just fine
  dplyr::select(YEAR_, ALARM_DATE, FIRE_NAME, geometry) %>% #select also seems to work and shouldn't mess up anything
  filter(YEAR_ == 2017) %>%  #so the problem is in 2017, but I can try to break it down into months with ALARM_DATE
  
  #so now if this valid works, we know the error is in the first half of the year
  #filter(ALARM_DATE > as.Date("2017-06-01")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #okay, this is all NA again

  #let's try the beginning of the year
  #filter(ALARM_DATE < as.Date("2017-06-01")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #THE BEGINNING IS FINE!!!

  #okay, we'll try october
  #filter(ALARM_DATE > as.Date("2017-10-01")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #After is all working!, so now the error must be between june 1st and october 1st

  #let's try september
  #filter(ALARM_DATE > as.Date("2017-09-01")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #september is fine

  #august?
  #filter(ALARM_DATE > as.Date("2017-08-01")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #august is fine
  
  #july?
  #filter(ALARM_DATE > as.Date("2017-07-01")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #broken again, so that means that the error is between july 1st and august 1st

  #just july
  #filter(ALARM_DATE > as.Date("2017-07-01") & ALARM_DATE < as.Date("2017-08-01")) 
  #still 148 fires

  #first half of July
  #filter(ALARM_DATE > as.Date("2017-07-01") & ALARM_DATE < as.Date("2017-07-15")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #the first half is clear

  #last 1/4 if July?
  #filter(ALARM_DATE > as.Date("2017-07-23") & ALARM_DATE < as.Date("2017-08-01")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #clear, so that means the error is between july 15th and july 23rd

  #try this window
  #filter(ALARM_DATE > as.Date("2017-07-14") & ALARM_DATE < as.Date("2017-07-24")) #%>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #error is in here, there are 31 observations, I will try manually checking
    #NOTHING obvious

  #filter(ALARM_DATE > as.Date("2017-07-16") & ALARM_DATE < as.Date("2017-07-19")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()

  #it's between 16 and 17
  #filter(ALARM_DATE > as.Date("2017-07-15") & ALARM_DATE < as.Date("2017-07-18")) #%>% 
  #st_is_valid() %>% 
  #as.data.frame()
  
#there's one :2017 2017-07-17 MURPHY  has multiple 3 point polys
  
  #filter(FIRE_NAME == "MURPHY" & ALARM_DATE > as.Date("2017-07-15") & ALARM_DATE < as.Date("2017-07-18")) %>% 
  #st_is_valid() %>% 
  #as.data.frame()
  #that's the one

  #now, let's try filtering it out
  filter(!(FIRE_NAME == "MURPHY" & ALARM_DATE > as.Date("2017-07-15") & ALARM_DATE < as.Date("2017-07-18"))) %>% 
  st_is_valid() %>% 
  as.data.frame()
## Reading layer `California_Wildland_Fire_Perimeters_All' from data source 
##   `/Users/Jason/github/fire_map_proj/data/California_Wildland_Fire_Perimeters_-All/California_Wildland_Fire_Perimeters_All.shp' 
##   using driver `ESRI Shapefile'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received an unexpected geometry. Either a polygon
## with interior rings, or a polygon with less than 4 points, or a non-Polygon
## geometry. Return arguments as a collection.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Geometry of polygon of fid 20145 cannot be translated to Simple
## Geometry. All polygons will be contained in a multipolygon.
## Simple feature collection with 21688 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13848330 ymin: 3833204 xmax: -12705610 ymax: 5255380
## Projected CRS: WGS 84 / Pseudo-Mercator
  #IT WORKS YEEEESSSSS, THIS TOOK AT LEAST 2 HOURS. I NEED A BREAK!
library(rmapshaper)
## Registered S3 method overwritten by 'geojsonlint':
##   method         from 
##   print.location dplyr
#fire map time
fire_data_path <- "../data/California_Wildland_Fire_Perimeters_-All/California_Wildland_Fire_Perimeters_All.shp"

#fire data has date format for dates
fire_data_perim <- st_read(fire_data_path) %>% 
  st_as_sf() %>% 
  #okay, so I' having issues with small polygons in the multipoly
  #I'll try unioning them
  dplyr::select(YEAR_, ALARM_DATE, FIRE_NAME, geometry) %>% 
  filter(YEAR_ > 2015)
## Reading layer `California_Wildland_Fire_Perimeters_All' from data source 
##   `/Users/Jason/github/fire_map_proj/data/California_Wildland_Fire_Perimeters_-All/California_Wildland_Fire_Perimeters_All.shp' 
##   using driver `ESRI Shapefile'
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, : GDAL
## Message 1: organizePolygons() received an unexpected geometry. Either a polygon
## with interior rings, or a polygon with less than 4 points, or a non-Polygon
## geometry. Return arguments as a collection.
## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## GDAL Message 1: Geometry of polygon of fid 20145 cannot be translated to Simple
## Geometry. All polygons will be contained in a multipolygon.
## Simple feature collection with 21688 features and 20 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -13848330 ymin: 3833204 xmax: -12705610 ymax: 5255380
## Projected CRS: WGS 84 / Pseudo-Mercator
#now select for fires in 2016
fire_data_2016 <- fire_data_perim %>% 
  filter(YEAR_ == 2016) %>% 
  #st_union() %>% #this eliminates the issue st_simplify was having with tiny polygons
  st_make_valid() %>% #valid polygons were triggering the between 0 and 4 points error when they clearly had more than 3 points
  rmapshaper::ms_simplify(keep = .1, keep_shapes = TRUE) %>% #ms_simplify does better to preserve topology and keeps 10% of the points
  st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet


#now select for fires in 2017
fire_data_2017 <- fire_data_perim %>% 
  filter(YEAR_ == 2017) %>%
  filter(!(FIRE_NAME == "MURPHY" & ALARM_DATE > as.Date("2017-07-15") & ALARM_DATE < as.Date("2017-07-18"))) %>% 
  st_make_valid() %>% 
  #st_union() %>% #this eliminates the issue st_simplify was having with tiny polygons
  rmapshaper::ms_simplify(keep = .1, keep_shapes = TRUE) %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet



#now select for fires in 2018
fire_data_2018 <- fire_data_perim %>% 
  filter(YEAR_ == 2018) %>% 
  st_make_valid() %>% 
  rmapshaper::ms_simplify(keep = .1, keep_shapes = TRUE) %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet

#now select for fires in 2019
fire_data_2019 <- fire_data_perim %>% 
  filter(YEAR_ == 2019) %>% 
  #might need to remove 'BELLA COLLINA GOLF COURSE SCL' and 'ORTEGA'
  filter(!(FIRE_NAME == 'ORTEGA' | FIRE_NAME == 'BELLA COLLINA GOLF COURSE SCL')) %>% 
  #yes, st_make_valid messed them up for some reason
  st_make_valid() %>% 
  rmapshaper::ms_simplify(keep = .1, keep_shapes = TRUE) %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet


#now select for fires in 2020
fire_data_2020 <- fire_data_perim %>% 
  filter(YEAR_ == 2020) %>% 
  st_make_valid() %>% 
  rmapshaper::ms_simplify(keep = .1, keep_shapes = TRUE) %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet

#now select for fires in 2021
fire_data_2021 <- fire_data_perim %>% 
  filter(YEAR_ == 2021) %>% 
  st_make_valid() %>% 
  rmapshaper::ms_simplify(keep = .1, keep_shapes = TRUE) %>% 
  st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet

#NO 2022 FIRE POLYGON DATA!!!!!

#now select for fires in 2022
#fire_data_2022 <- fire_data_perim %>% 
  #filter(YEAR_ == 2022) %>% 
  #st_union() %>% #this eliminates the issue st_simplify was having with tiny polygons
  #t_make_valid() %>% 
  #st_simplify(dTolerance = 200) %>% 
  #st_transform(crs = "+proj=longlat +datum=WGS84") #needed for projecting on leaflet

#create ggplot
ggplot()+
  #color the objects
  geom_sf(data = fire_data_2016, fill = "#454742", col = "#8d8c8b", size = .5)+
  #create the theme for the background. I'm reusing an old theme
  theme(panel.background = element_rect(fill = "#272727",
                                        colour = "#272727",
                                        size = .5,
                                        linetype = "dashed"),
        panel.grid.major = element_line(size = .5,
                                        color = "#3d3d3d",
                                        linetype = "solid"))+
  #title and subtitle. I might need another legend position for credits.
  theme(legend.position = "hide")+
    labs(title = "TITLE",
         subtitle = "SUBTITLE",
         x = "",
         y ="")

Now that the fire data is tidy and ready for the map, I’ll tidy up the precip data

#do I need raster libraries
library(raster)
library(rasterVis)
## Loading required package: terra
## terra version 1.3.4
## 
## Attaching package: 'terra'
## The following object is masked from 'package:rgdal':
## 
##     project
## The following object is masked from 'package:dplyr':
## 
##     src
## Loading required package: lattice
## Loading required package: latticeExtra
## 
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
## 
##     layer
library(ggplot2)

#first get the us boundary CA state poly
ca_state_poly <- us_boundaries(states = "California")
## old-style crs object detected; please recreate object with a recent sf::st_crs()
## old-style crs object detected; please recreate object with a recent sf::st_crs()
#now get the precip data
GDALinfo("../data/nws_precip_ytd_20220101_conus.tif")
## Warning: GDAL support is provided by the sf and terra packages among others
## Warning: GDAL support is provided by the sf and terra packages among others
## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver

## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver

## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver

## Warning in GDALinfo("../data/nws_precip_ytd_20220101_conus.tif"): statistics not
## supported by this driver
## Warning: GDAL support is provided by the sf and terra packages among others
## rows        881 
## columns     1121 
## bands       4 
## lower left origin.x        -1904912 
## lower left origin.y        -7619987 
## res.x       4763 
## res.y       4763 
## ysign       -1 
## oblique.x   0 
## oblique.y   0 
## driver      GTiff 
## projection  +proj=stere +lat_0=90 +lat_ts=60 +lon_0=-105 +x_0=0 +y_0=0 +R=6371200
## +units=m +no_defs +type=crs 
## file        ../data/nws_precip_ytd_20220101_conus.tif 
## apparent band summary:
##    GDType hasNoDataValue   NoDataValue blockSize1 blockSize2
## 1 Float32           TRUE -3.402823e+38        256        256
## 2 Float32           TRUE -3.402823e+38        256        256
## 3 Float32           TRUE -3.402823e+38        256        256
## 4 Float32           TRUE -3.402823e+38        256        256
## apparent band statistics:
##          Bmin       Bmax Bmean Bsd
## 1 -4294967295 4294967295    NA  NA
## 2 -4294967295 4294967295    NA  NA
## 3 -4294967295 4294967295    NA  NA
## 4 -4294967295 4294967295    NA  NA
## Metadata:
## AREA_OR_POINT=Area 
## creation_time=20220112T010818 
## data_time=20220101
precipitation_conus_path <- "../data/nws_precip_ytd_20220101_conus.tif"

#yearly precip paths
#the date on the tif file is the 1st of the next year, when the data was presumably compiled
#2016
precip_2016_path <- "../data/raw_precip_data/nws_precip_ytd_20170101_geotiff_year_2016/nws_precip_ytd_20170101_conus.tif"

#2017
precip_2017_path <- "../data/raw_precip_data/nws_precip_ytd_20180101_geotiff_year_2017/nws_precip_ytd_20180101_conus.tif"

#2018
precip_2018_path <- "../data/raw_precip_data/nws_precip_ytd_20190101_geotiff_year_2018/nws_precip_ytd_20190101_conus.tif"

#2019
precip_2019_path <- "../data/raw_precip_data/nws_precip_ytd_20200101_geotiff_year_2019/nws_precip_ytd_20200101_conus.tif"

#2020
precip_2020_path <- "../data/raw_precip_data/nws_precip_ytd_20210101_geotiff_year_2020/nws_precip_ytd_20210101_conus.tif"

#2021
precip_2021_path <- "../data/raw_precip_data/nws_precip_ytd_20220101_geotiff_year_2021/nws_precip_ytd_20220101_conus.tif"


#create the raster 2022 rainfall
precip_conus <- raster(precipitation_conus_path, band = 2) %>% 
  projectRaster(crs = 4326) %>% 
  #crop california data
  #crop(ca_state_poly) %>% 
  #select only cali borders
  raster::intersect(ca_state_poly) %>% 
  #mask the state polygon
  raster::mask(ca_state_poly) %>% 
  #convert to pts
  rasterToPoints(spatial = TRUE) %>% 
  #then to df
  data.frame()
## Warning: PROJ support is provided by the sf and terra packages among others
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
#create the raster for leaflet
precip_conus_leaf <- raster(precipitation_conus_path, band = 2) %>% 
  projectRaster(crs = 4326) %>% 
  #crop california data
  #crop(ca_state_poly) %>% 
  #select only cali borders
  raster::intersect(ca_state_poly) %>% 
  #mask the state polygon
  raster::mask(ca_state_poly)
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
#rasters by year for the leaflet map 
precip_2016 <- raster(precip_2016_path, band = 2) %>% #band 2 has the precip data in inches
  projectRaster(crs = 4326) %>% #wgs 84 prjections
  raster::intersect(ca_state_poly) %>% #select only cali state extent
  raster::mask(ca_state_poly) #mask only the california data
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
precip_2017 <- raster(precip_2017_path, band = 2) %>% #band 2 has the precip data in inches
  projectRaster(crs = 4326) %>% #wgs 84 prjections
  raster::intersect(ca_state_poly) %>% #select only cali state extent
  raster::mask(ca_state_poly) #mask only the california data
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
precip_2018 <- raster(precip_2018_path, band = 2) %>% #band 2 has the precip data in inches
  projectRaster(crs = 4326) %>% #wgs 84 prjections
  raster::intersect(ca_state_poly) %>% #select only cali state extent
  raster::mask(ca_state_poly) #mask only the california data
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
precip_2019 <- raster(precip_2019_path, band = 2) %>% #band 2 has the precip data in inches
  projectRaster(crs = 4326) %>% #wgs 84 prjections
  raster::intersect(ca_state_poly) %>% #select only cali state extent
  raster::mask(ca_state_poly) #mask only the california data
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
precip_2020 <- raster(precip_2020_path, band = 2) %>% #band 2 has the precip data in inches
  projectRaster(crs = 4326) %>% #wgs 84 prjections
  raster::intersect(ca_state_poly) %>% #select only cali state extent
  raster::mask(ca_state_poly) #mask only the california data
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
precip_2021 <- raster(precip_2021_path, band = 2) %>% #band 2 has the precip data in inches
  projectRaster(crs = 4326) %>% #wgs 84 prjections
  raster::intersect(ca_state_poly) %>% #select only cali state extent
  raster::mask(ca_state_poly) #mask only the california data
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
ggplot()+
  #color the objects
  geom_sf(data = ca_state_poly, fill = "#454742", col = "#8d8c8b", size = .5)+
  #create the theme for the background. I'm reusing an old theme
  theme(panel.background = element_rect(fill = "#272727",
                                        colour = "#272727",
                                        size = .5,
                                        linetype = "dashed"),
        panel.grid.major = element_line(size = .5,
                                        color = "#3d3d3d",
                                        linetype = "solid"))+
  #title and subtitle. I might need another legend position for credits.
  geom_raster(data = precip_conus, aes(x=x, y=y, fill = nws_precip_ytd_20220101_conus))+
  scale_fill_gradientn(colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"))+
  theme(legend.position = "hide")+
    labs(title = "TITLE",
         subtitle = "SUBTITLE",
         x = "",
         y ="")

FINALLY!!!! I had to use the raster::mask fuction to do what I wanted cause raster::intersect only selected by the extent, either way, I have the tidy raster data I was looking for. The next step should be setting up the leaflet map. I think the best way to go about this will be to create the leaflet map I want in this RMD, but then create a new RMD that will be cleaner and just show the leaflet map and have this rmd in a different folder showing my process.

I’m going to make a chunk to work on the graphs here

so the graph I want is going to be a graph of the calculated yearly number of people who are ‘affected’ by fires by selecting counties with large fires in them or many small fires. I’m gonna need to start by just simply selecting fires based on size and figure out the smaller fires issue later.

#select fires based on size

#might need specific county
cali_counties_sonoma <- cali_counties_final %>% 
  filter(county_name == "Sonoma") %>% 
  st_make_valid()

#let's see what the raw fire data has that could help me
fire_data_select_test <- fire_data_2017 %>% #okay, from here I'm going to try to add a column containing the county names the fire intersects
  st_transform(crs = "+proj=longlat +datum=WGS84") %>%
  st_make_valid() %>% #had to make the fire data clean because of the conversion from geos to s2 and duplicate points
  st_intersection(cali_counties_sonoma)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
ggplot()+
  #color the cali counties and their borders
  geom_sf(data = cali_counties_sonoma, fill = "#454742", col = "#8d8c8b", size = .5)+
  #create the theme for the background. I'm reusing an old theme
  theme(panel.background = element_rect(fill = "#272727",
                                        colour = "#272727",
                                        size = .5,
                                        linetype = "dashed"),
        panel.grid.major = element_line(size = .5,
                                        color = "#3d3d3d",
                                        linetype = "solid"))+
  #title and subtitle. I might need another legend position for credits.
  geom_sf(data = fire_data_select_test, fill= "red", col = "red")#check the fire data

  theme(legend.position = "hide")+
    labs(title = "TITLE",
         subtitle = "SUBTITLE",
         x = "",
         y ="")
## List of 5
##  $ legend.position: chr "hide"
##  $ x              : chr ""
##  $ y              : chr ""
##  $ title          : chr "TITLE"
##  $ subtitle       : chr "SUBTITLE"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

Okay, so that works with just a singular county, but does it work with all counties?

#select fires based on size

#might need specific county
cali_counties_valid <- cali_counties_final %>% 
  st_make_valid()

#need to strip the geometry from counties for joining
counties_mini_test <- cali_counties_valid %>% 
  mutate(county_area = st_area(cali_counties_valid$geometry)) %>% 
  as.data.frame() %>% 
  dplyr::select(county_name, county_population, county_area)

#let's see what the raw fire data has that could help me
fire_data_select_test_full <- fire_data_2020 %>% #okay, from here I'm going to try to add a column containing the county names the fire intersects
  st_transform(crs = "+proj=longlat +datum=WGS84") %>%
  st_make_valid() %>% #had to make the fire data clean because of the conversion from geos to s2 and duplicate points
  st_intersection(cali_counties_valid) %>% #okay, so this gave me what I wanted, it attached the county name and population count to the fires
  group_by(county_name) %>% #group the fires by county name so I can operate on them
  summarize(geometry = st_union(geometry)) %>% #summarize is being weird, but I can just add the county pop after
  left_join(counties_mini_test, by = c('county_name' = 'county_name')) %>%  #okay, hopefully, geometry not being the last column isn't a problem
  # 2020 data error: Loop 5 is not valid: Edge 14 is degenerate (duplicate vertex)
  st_make_valid() %>% #have to call again because the prior steps after make_valid can mess stuff up

  mutate(fire_area = st_area(geometry)) %>%  #add fire area
  mutate(fire_coverage = as.numeric(fire_area/county_area)*100) %>%  #percentage of county area covered by fires
  mutate(fire_impacted = ifelse(fire_coverage > 1, "impacted", "not impacted")) %>% 
  group_by(fire_impacted) %>% 
  summarize(impacted_pop = sum(as.numeric(county_population))) #now I have the number of impacted people in the impacted_pop column
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
#I think I'm going to judge that any county where over 1% of the county area was covered by fires contitutes the county being
#'impacted by wildfires' for the purposes of this qualitiative analysis. I could also say that about 8% is 'seriously impacted'
#since that's about the fire coverage for sonoma county in 2017.

#okay, I've got a column of impacted fires, next step will be figuring out how I want to put it in the leaflet map.

#the next step is to create a function to automate the impact plot thing and create a data frame with
#year on the x-axis and impacted_pop on the y-axis for the ggplot.

   
  

plot(fire_data_select_test_full)

test_graph <- ggplot()+
  #color the cali counties and their borders
  geom_sf(data = cali_counties_valid, fill = "#454742", col = "#8d8c8b", size = .5)+
  #create the theme for the background. I'm reusing an old theme
  theme(panel.background = element_rect(fill = "#272727",
                                        colour = "#272727",
                                        size = .5,
                                        linetype = "dashed"),
        panel.grid.major = element_line(size = .5,
                                        color = "#3d3d3d",
                                        linetype = "solid"))+
  #title and subtitle. I might need another legend position for credits.
  geom_sf(data = fire_data_select_test_full, fill= "red", col = "red")#check the fire data
  theme(legend.position = "hide")+
    labs(title = "TITLE",
         subtitle = "SUBTITLE",
         x = "",
         y ="")
## List of 5
##  $ legend.position: chr "hide"
##  $ x              : chr ""
##  $ y              : chr ""
##  $ title          : chr "TITLE"
##  $ subtitle       : chr "SUBTITLE"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi FALSE
##  - attr(*, "validate")= logi TRUE

Before I make the graph in ggplot, I need to assemble the dataframe that holds the year and pop data

#create a function to find the impacted pop
find_impacted_pop <- function(fire_data_func){
  
  impacted_pop_func <- fire_data_func %>% 
    st_transform(crs = "+proj=longlat +datum=WGS84") %>%
  
    st_make_valid() %>% #had to make the fire data clean because of the conversion from geos to s2 and duplicate points
    st_intersection(cali_counties_valid) %>% #okay, so this gave me what I wanted, it attached the county name and population count to the fires
    group_by(county_name) %>% #group the fires by county name so I can operate on them
    summarize(geometry = st_union(geometry)) %>% #summarize is being weird, but I can just add the county pop after
    left_join(counties_mini_test, by = c('county_name' = 'county_name')) %>%  #okay, hopefully, geometry not being the last column isn't a problem
    st_make_valid() %>% #call again to clean up before st_area
    mutate(fire_area = st_area(geometry)) %>%  #add fire area
    mutate(fire_coverage = as.numeric(fire_area/county_area)*100) %>%  #percentage of county area covered by fires
    mutate(fire_impacted = ifelse(fire_coverage > 1, "impacted", "not impacted")) %>% 
    group_by(fire_impacted) %>% 
    summarize(impacted_pop = sum(as.numeric(county_population)))
    
  return(impacted_pop_func)
}

#run each year through the function
fire_impacted_2016 <- find_impacted_pop(fire_data_2016)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
fire_impacted_2017 <- find_impacted_pop(fire_data_2017)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
fire_impacted_2018 <- find_impacted_pop(fire_data_2018)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
fire_impacted_2019 <- find_impacted_pop(fire_data_2019)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
fire_impacted_2020 <- find_impacted_pop(fire_data_2020) #problem with this one
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
fire_impacted_2021 <- find_impacted_pop(fire_data_2021)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
#assemble the data frame with year and pop impacted
impacted_population_df <- data.frame(
  YEAR = c("2016","2017","2018","2019", "2020", "2021"),
  IMPACTED = c(fire_impacted_2016$impacted_pop[1],
               fire_impacted_2017$impacted_pop[1],
               fire_impacted_2018$impacted_pop[1],
               fire_impacted_2019$impacted_pop[1],
               fire_impacted_2020$impacted_pop[1],
               fire_impacted_2021$impacted_pop[1])
)

Now I’m going to work on the graph itself.

library(scales)
## 
## Attaching package: 'scales'
## The following object is masked from 'package:terra':
## 
##     rescale
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
options(scipen = 999)
  
impact_graph <- ggplot(data = impacted_population_df, aes(x = YEAR, y = IMPACTED))+
  geom_col(col = "#ff6400", fill = "#ff6400")+
  labs(title = "People Impacted by Wildfires",
       x = "Year",
       y = "Impacted Population",
       subtitle = "Counties with over 1% fire area coverage considered impacted.") +
  scale_y_continuous(labels = scales::label_number_si())+
  theme_minimal()+
  theme(axis.text = element_text(face="bold", size = rel(2)),
        axis.title = element_text(face="bold", size = rel(1)),
        title = element_text(face="bold", size = rel(2)),
        axis.ticks = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.major.y = element_line(color = "black"),
        panel.grid.minor.y = element_line(color = "grey"),
        plot.subtitle = element_text(size = rel(.55))
        )
  

impact_graph

#time to try the leaflet map
library(leaflet)
library(leaflet.extras)
library(leafpop)


#html code and stuff to display the county info
county_popup_text <- paste("<div class='leaflet-popup-scrolled' style='max-width:600px;max-height:300px'><b>",
                           
                           '  <h2 align="center"><strong>', cali_counties_final$county_name, "</strong></h2></a><br>",
               "</b>",
               
               ' <b>County Population: </b> ',
                           
                            cali_counties_final$county_population,'<br><br><center><img src="')


graph_icon <- makeIcon(iconUrl = "../data/bar-chart-outline.svg", 
                       iconWidth = 50,
                       iconHeight = 50)




map = leaflet() %>%
  addTiles(group = "OSM") %>% #open street map baselayer
  addProviderTiles(providers$CartoDB.Positron, group = "Carto Light") %>% #cartodb light basemap
  addProviderTiles(providers$CartoDB.DarkMatter, group = "Carto Dark") %>%  #carto db dark basemap
  addPolygons(data = cali_counties_final, 
              fillOpacity = 0, 
              weight = 1, 
              color = "grey", 
              group = "County Borders",
              popup = county_popup_text) %>% #county border polys
  

  
  addPolygons(data = ca_places_bound, 
              color = "grey", 
              fillOpacity = .4, 
              weight = 0, 
              group = "Pop Centers",
              popup = ca_places_bound$NAMELSAD) %>% 
  
  #2016 fires
  addRasterImage(x = precip_2016, 
                 opacity = .5, 
                 colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"), 
                 group = "Fires in 2016") %>% 
  addPolygons(data = fire_data_2016, 
              opacity = .5, 
              color = "red", 
              group = "Fires in 2016",
              popup = fire_data_2016$FIRE_NAME) %>% 
  
  #2017 fires
  addRasterImage(x = precip_2017, 
                 opacity = .5, 
                 colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"), 
                 group = "Fires in 2017") %>% 
  addPolygons(data = fire_data_2017, 
              opacity = .5, 
              color = "red", 
              group = "Fires in 2017",
              popup = fire_data_2017$FIRE_NAME) %>% 
  
  #2018 fires
  addRasterImage(x = precip_2018, 
                 opacity = .5, 
                 colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"), 
                 group = "Fires in 2018") %>% 
  addPolygons(data = fire_data_2018, 
              opacity = .5, 
              color = "red", 
              group = "Fires in 2018",
              popup = fire_data_2018$FIRE_NAME) %>% 
  
  #2019 fires
  addRasterImage(x = precip_2019, 
                 opacity = .5, 
                 colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"), 
                 group = "Fires in 2019") %>% 
  addPolygons(data = fire_data_2019, 
              opacity = .5, 
              color = "red", 
              group = "Fires in 2019",
              popup = fire_data_2019$FIRE_NAME) %>% 
  
  #2020 fires
  addRasterImage(x = precip_2020, 
                 opacity = .5, 
                 colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"), 
                 group = "Fires in 2020") %>% 
  addPolygons(data = fire_data_2020, 
              opacity = .5, 
              color = "red", 
              group = "Fires in 2020",
              popup = fire_data_2020$FIRE_NAME) %>% 
  
  #2021 fires
  addRasterImage(x = precip_2021, 
                 opacity = .5, 
                 colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"), 
                 group = "Fires in 2021") %>% 
  addPolygons(data = fire_data_2021, 
              opacity = .5, 
              color = "red", 
              group = "Fires in 2021",
              popup = fire_data_2021$FIRE_NAME) %>% 

  #layer controls
  addLayersControl(
    baseGroups = c("OSM", "Carto Light", "Carto Dark"),
    overlayGroups = c("County Borders", 
                      "Precipitation",
                      "Pop Centers",
                      "Fires in 2016",
                      "Fires in 2017",
                      "Fires in 2018",
                      "Fires in 2019",
                      "Fires in 2020",
                      "Fires in 2021"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
  #create legend
  addLegend(position = "bottomleft", 
            colors = c("#fff8bf", "#6e6a4b",  "#141b57", "#000724"), 
            labels = c("Low Rainfall", "", "", "High Rainfall"), 
            group = "Legend") %>% 
  addSearchOSM() %>% 
  addFullscreenControl() %>% 
  hideGroup("Fires in 2016") %>% 
  hideGroup("Fires in 2017") %>% 
  hideGroup("Fires in 2018") %>% 
  hideGroup("Fires in 2019") %>% 
  hideGroup("Fires in 2020") %>% 
  addMiniMap() %>% 
  addMarkers(lng=-117, lat=39, popup = popupGraph(impact_graph, width = 450, height = 300), icon = graph_icon)
## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others

## Warning: PROJ support is provided by the sf and terra packages among others
  #implementation: https://rdrr.io/github/r-spatial/leafpop/man/addPopupGraphs.html
  #I still need to make the actual data for the graph first
  #addPopupGraphs()
  #I decided to just use a ggplot and put it in a popup
  
  

map

Okay, great, this map is looking pretty mediocre, I can worry more about the aesthetics later, for now I have the most base layer finished. Now I have to figure out how to make an interactive switching map.

Right, so now I’m just remembering that I only processed the precipitation raster for 2022, so now I need to go back up and finish processing all the year layers.

Alrighty, so the map itself is done, now I just need to make the graphs.

Welp, it’s done now. I have more ideas of what I could add to the map, but at this point, I think it’s about time to move on to other things. I could add different colors to the counties based on whether the county is ‘impacted’ but that’s not really necessary and this map is already detailed enough. I’m happy with the outcome and if I were given this map as a task to complete it would certainly be ready for review by the client so I’m going to leave it in the state it’s in now.

The next step is to duplicate this rmd and remove all of the unnecessary bits of code and commentary so that I can turn it into a nice webpage for my website. This rmd will likely be linked as a ‘journal’ of sorts, documenting my creation process as I built this map from scratch.